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Computer methods have recently been applied to a large class 
of steady-state collisionless plasma sheath problems, including, for 

example, moving satellites in the ionosphere,^ stationary plasma 

v ( 6 “ 8 ) • • (9) , (10) 

probes, ion engines, and plasma diodes. These problems 


involve numerical solutions of the Poisson equation. One aspect of 
such solutions which can cause difficulties is that in many cases ^ 


one boundary is at infinity, where the potential and net charge density 
vanish. Because of the limitations of digital computers, the potential 

and density descriptions are restricted to a finite portion of space, so 
that the boundary condition corresponding to "infinity" must be simulated 
by a necessarily artificial finite condition. Another aspect is that 
the density is generally a non-linear functional of the potential 
distribution, and that, with specified boundary conditions, an iteration 
procedure must be used in order to obtain a self-consistent solution. 

The stability of the iteration procedure depends, not only on the 
particular algorithm, but also on the artificial boundary condition 
employed, as will be shown later. Even in finite problems bounded by 




electrodes,^ the stability of the iteration is related to the 

magnitude of the electrode separation.^ 

This paper is concerned with the numerical effects of employing 
certain artificial finite boundary conditions and iterative procedures 
for self-consistent Poisson problems. Studies have been made with 
spherical probe and plasma diode models, with the goal of developing 
economical procedures for more general problems. Results obtained by 
the use of such models are not expected to depend strongly on their 
one-dimens iona 1 i ty . A third aspect of the plasma sheath problem, which 
is numerically non-trivial but will not be discussed, is that of finding 
the density when the potential distribution is given, in the presence of 
absorbing boundaries and non-isotropic velocity distributions. Alter- 
native approaches for this calculation in a general problem are suggested 
in Refs. 2-6 and 12. 

A common approach to the artificial finite boundary has been 
simply to set the potential to zero on a boundary surface as far "out" 
as possible, within computer limitations. "Float ing-potential " 

boundary conditions have also been used, ^ ^ in which linear relation- 

ships have been assumed between the potential and its gradient. The 
question of how well a finite boundary condition approximates the "true" 

( 3 - 8 ) 

infinite boundary condition has been treated recently. In these 

invest iga tions , sequences of Poisson problems with identical boundary 

conditions were solved, where the boundary was successively moved 

outward until no further changes occurred, either in the potential 

(3-4) 

distribution in the vicinity of the satellite, or in the collected 


/ 


2 



( 5 - 8 ) 

probe current'.. The solution obtained in this way may be taken to 

represent that for the "infinite 11 boundary condition, and will be said 

to have a "stationary" property with respect to the boundary position. 

Computational costs are proport Iona 1 to the size of the region 

enclosed by the boundary. If the boundary can be moved inward as a 

result of changing the boundary condition, without changing the stationary 

property of the solution, the calculation becomes more economical. In 

the case of a stationary spherical probe, for example, the potential is 

known to obey an inverse-square asymptotic lav;, and computer procedures 

( 7 - 8 ) 

employing floating boundary conditions based on this law v have been 

found to be far more economical than procedures in which the potential 

( 4 ) 

is set to zero. In the case of a moving satellite or a circular planar 

probe, ^ ^ where axially symmetric potential distributions occur and no 

theoretical asymptotic analysis is available, various ad ho c floating- 

potential relationships, based on intuition, have similarly been found 

( 4 ) 

to be more economical than setting the potential to zero. Taylor, 

for example, has employed an exponential law as well as certain other 
special relationships, while Parker and Whipple^ ^ have employed a 
dipole law. In view of these results, it seems worthwhile to look for 
generally efficient assumptions. For example, for the spherical probe 

we have found that setting the gradient of the potential to zero at the 

( 13 ) 

boundary is as economical as using the asymptotically correct inverse 

square law. Calculations will be presented which will compare floating 
conditions with the zero-potential condition. 
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The iterative phase of the procedure for obtaining self-consistent 
numerical solutions may be discussed in terms of the Poisson matrix equa- 
tion in dimensionless form, namely: 


L <t> = - p(^) (1) 

In (1), L is the Laplacian matrix operator resulting from the differencing 
of the Laplacian on a grid in conf igura tiona 1 space. The components of 
the vectors i and p are the values of the potential and the net charge 
density, respectively, at the grid points. A general iterative procedure 
may be based on the mixing or coupling of density iterates, namely: 


l i 


n-t-l 
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m-0 


where 0^ denotes the n-th iterate for Alternatively, of course, one 

may instead couple the potential iterates. The coefficients b^ determine 

the stability of the iteration.* 

Using the spherical probe model, we find that the iteration 

converges when the boundary radius R is less than a certain critical 

value , and diverges when R is greater than . Let denote the 

boundary radius such that the probe current is stationary when R > R^ . 

The iteration will converge to the desired solution (with respect to 

the current) when R > R . With "direct" (uncoupled) iteration, in 

c s 

which b = 8 (Kronecker delta), R is usually less than R . Now R is 
m mn > c s s 

fixed by the chosen boundary condition, but R^ can be made to become 



In the iteration scheme with mixing parameter CX < 1, defined by 


n-bl 


W n+1 - F_, - - op CfJ + (1 - a) F„ (3) 


n 


(where F > =0 so that 0^ is the Laplace solution), the stability of the 

iteration (3) is found to increase as a decreases. For any value of R^, 

one can apparently force convergence, i.e., cause R to exceed R^ , if one 

chooses a sufficiently small value of a. Equation (3) implies that b 

in (2) is given by a(l - a) n m , so that for small a the coupling is 

approximately equally distributed over a large number of past iterates. 

The use of such an iteration corresponds to replacing the problem by an 

equivalent time-dependent parabolic heat-diffusion problem in difference 

form, in which the time index is analogous to the iteration index, and 

( 2 ) 

the steady-state solution is the desired one. The use of a small 

value for a corresponds to the introduction of high damping in the 
heat-diffusion problem. The required number of iterations, and therefore 
computer time, is roughly inversely proportional to a. 

Various versions of the iteration scheme (3) have been 
usec!.^ ** ® Laframboise ^ has found that significant gains 

can result from additional sophistication, such as the use of an empirical 
diagonal matrix for a rather than a scalar, that is, a different mixing 
parameter for each grid point. We are pursuing this promising approach 
and hope to report new results. 

One-dimensional plasma diode problems also afford convenient 
models for studying iterative stability. ^ ^ Numerical instability 
seems to be connected with the existence of an extended region of space 
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where the charge density is small, in the vicinity of the point where it 

changes sign. As the electrode spacing increases, the extent of the regime 

of small charge density grows, and the parameter a must be reduced for 

convergence. If, in the plasma diode problem, x^e restrict attention to 

a region in the vicinity of the point where the charge density changes 

sign, and if x*/e assume a linear relation between p and <p in this region, 

then it can be shown analytically that CC must be smaller than a number 
2 2 

proportional to X /I) , where X is the Debye length, and D is the extent 

of the region of linearity. If D is taken to be the electrode separation 

in an empirical formula, the data of Prince and Jeffries^^ indicates 

a constant of proportionality between 1 and 10 for large values of D/X. 

It is interesting to note that in satellite problems with high ion Mach 
(2-3) 

numbers, iterative stability is found to increase with increasing 

Mach number. One can probably define in this case an effective Debye 

length, which is based on the ion energy and is therefore large. 

In our experiments with the spherical probe model we have 

found that, for interesting values of the physical parameters, the scheme 

of Eq. (3) requires a to be of the order of 1/100. This small value 

implies that many iterations v. 7 ill be required, i.e., of the order of 
(14) 

hundreds. Details of these experiments will be presented* 

Based on the parabolic nature of the iteration problem, it 
should be possible to devise more sophisticated iteration schemes which 
will reduce the number of iterations required and thereby reduce computer 
t ime . 
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